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In current practice, it is often difficult to draw firm conclusions about turbulence model ac- 
curacy when performing multi-code CFD studies ostensibly using the same model because 
of inconsistencies in model formulation or implementation in different codes. This paper de- 
scribes an effort to improve the consistency, verification, and validation of turbulence mod- 
els within the aerospace community through a website database of verification and validation 
cases. Some of the variants of two widely-used turbulence models are described, and two inde- 
pendent computer codes (one structured and one unstructured) are used in conjunction with 
two specific versions of these models to demonstrate consistency with grid refinement for sev- 
eral representative problems. Naming conventions, implementation consistency, and thorough 
grid resolution studies are key factors necessary for success. 


1 . Introduction 

Computational fluid dynamics (CFD) codes that solve the Reynolds-averaged Navier-Stokes (RANS) 
equations are applied regularly at government laboratories, universities, and companies throughout the 
world to predict flow fields in and around complex configurations. A major component of most RANS 
applications is the turbulence model, which attempts to model the effects of turbulent fluctuations on 
the mean flow. Many types of turbulence models exist, ranging from algebraic models through linear 
one- and two-equation and all the way to nonlinear seven-equation full Reynolds-stress. One thing that 
all models have in common is that they are imperfect: they are unable to reproduce all features of all 
classes of flows. Unfortunately, it is often difficult to isolate whether the cause of any particular failure of 
a given code to match an experiment is due to the turbulence model itself, or some other reason (e.g., 
numerical errors, boundary conditions, geometry). 

One disturbing issue that has surfaced over the years is the fact that different codes that use osten- 
sibly the same turbulence model often do not approach the same result for the same case as the grid is 
refined (see, e.g., Vassberg et al.[1]). Beside inadequately refined grids, two of the contributing causes 
to this state of affairs is that different codes either knowingly or unknowingly use different versions of 
a turbulence model, or that a particular implementation has not been thoroughly verified. It can some- 
times be difficult to implement turbulence models from the open literature because of poor reporting 
practices [2]. Verification of turbulence models can also be problematic. One way is through the method 
of manufactured solutions (MMS), but only limited solutions are available for turbulence (see, e.g., Eca 
et al.[3]). Demonstrating consistency using the same model in different codes can provide some level 
of assurance [4], but optimally the testing should be conducted by independent groups [5] to minimize 
the possibility of repeating the same mistakes. 

In an effort to help improve the consistency, verification, and validation of turbulence models within 
the aerospace community, NASA has recently established a website to provide a central location where 
widely-used RANS turbulence models are described, grids and cases are provided, and detailed results 
are given (http://turbmodels.larc.nasa.gov, cited 19 May 2009). As it grows, this site should serve as 
a resource to the CFD community for verifying turbulence model implementations. This capability is 
made possible by providing simple test cases and grids, along with sample results (including complete 
grid convergence studies) from one or more previously-verified codes. Furthermore, by listing accepted 
versions of the turbulence models as well as published variants, this site establishes naming conventions 
in order to help avoid confusion when comparing results from different codes. 

The above-mentioned website is also a part of a larger effort, involving a discussion group of the 
American Institute of Aeronautics and Astronautics (AIAA), to establish a database and standards for 
turbulence modeling assessment. A similar goal was initially outlined at a turbulence modeling workshop 
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held in 2001 , as described in Rubinstein et al.[6]. The current paper primarily describes the verification 
part of the effort, and stresses the importance of consistency in implementation as well as conducting 
thorough grid convergence studies. 


2. Examples of Turbulence Model Variants 

Even among the widely-used one-equation Spalart-Allmaras (SA) [7] and two-equation Menter shear 
stress transport (SST) [8] k-w models, several variants exits. Sometimes, these variations are glossed 
over or ignored when reporting results. This can result in misunderstanding and/or misinterpretation 
of trends. The current effort seeks to fully document and label many turbulence model variants, thus 
establishing a basis from which future identification and cross-correlation of models and their results will 
be possible. This basis is a key component in the effort to achieve consistency between different codes. 

The one-equation SA model is written in terms of the turbulence quantity v\ 
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where a description of each of the terms is not given here, but can be found in the original reference [7]. 
The quantity v is related to the eddy viscosity by: 
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where v is the molecular kinematic viscosity and c vl = 7.1. The so-called SA-la version of this model in- 
cluded an additional trip term, but by common practice most codes do not include it. The recommended 
far field boundary condition [9, 10] for Eq. (1) is in the range of 3^ < v far field, < 5i/oo, and at the body 
Vwaii = 0. For all results reported here, v far field, = 

Some programmers leave out the / t2 term in Eq. (1), since it was originally associated with the use 
of the trip term. This version of the model, termed here SA-noft2, can be found for example in Eca et 
al.[3]. Other versions of the Spalart-Allmaras model include SA-RC [11], SA-Catris [12], SA-Edwards 
[13], and SA-salsa [14], 

The two-equation SST model is written in terms of the two turbulence quantities k and u>. The form 
is: 
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The eddy viscosity is given by: 


d\k 

At = P 7 (5) 
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where ai = 0.31, Ct is the magnitude of vorticity, and F 2 is a blending function. The production term V is 
given by 
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The recommended wall boundary conditions are: k wa ii = 0, u wa ii = 60i//[/3i(Adi) 2 ]. A complete 
description of each of the terms in the standard SST equations can be found in Menter [8]. 

A common variant of the SST model is termed SST-V, for which the production term makes use of 
the local magnitude of vorticity Q: 


v = mti 1 - \p k8 a^: ( 9 ) 

This vorticity source term is often close to the exact source term in boundary layer flows [15], and its 
use can avoid some numerical difficulties sometimes associated with the use of the exact source term. 
Other variants of SST include SST-2003 [16] and SST-sust [10]. The main change in SST-2003 is a 
modification to the definition of y t , Eq. (5): 


P max(a 1 aj,SF 2 ) 1 } 

where S = y / 25' ij -5'y. The SST-sust model is identical to the original SST model except for the addition 
of constant sustaining terms /3*w /or/ieW fc /or/ieW and /3w} arfield in Eqs. (3) and (4), respectively. In 
the free stream, these have the effect of exactly canceling the destruction terms if k = kf ar field and 
uj = ujfarfieid, thus preserving free stream levels of k and u>. Inside the boundary layer, the terms are 
generally orders of magnitude smaller than the destruction terms for reasonable free stream turbulence 
levels (say, Tu = 1 % or less), and therefore have little effect. The main advantage of the use of sustaining 
terms is that they remove ambiguity associated with far field boundary conditions. For SST (as well as 
other similar two-equation models), there is free stream decay of k and u, and local ambient levels 
can vary depending on distance to the outer boundary. For SST-sust, k and u do not decay in the 
free stream. The recommended far field boundary conditions for SST-sust are: kf ar f ie i d = 10 ~ 6 U^ 
and ujfarfieid = SUoo/L. The far field boundary conditions for SST, SST-V, and SST-2003 tend to vary 
depending on the code and application. For all results reported here, k far f ie id = 9 x 10 _9 a^ and 

kfarfieid = 1 X 10 Poo&oo/ Poo‘ 

At the time of this writing, only the SA and SST models and their variants have been posted to the 
website. As the website database grows, it is planned to eventually include other models, particularly 
those that see wide use in the aerospace community. 


3. Results 

Sample results are first described for flow over a zero-pressure-gradient (ZPG) flat plate using two 
different compressible Navier-Stokes CFD codes, CFL3D [17] and FUN3D [18]. The former is a cell- 
centered structured code and the latter node-centered unstructured. These two codes have both un- 
dergone some level of verification using MMS for turbulent flows [4], The flow conditions are M = 0.2, 
Re = 10 million based on the full length of the plate ( L = 2), or Re = 5 million per unit length. Although 
not shown, when using SA and SST-V both codes yield reasonable results in terms of reproducing well- 
known ZPG law-of-the-wall behavior and wall skin friction for which correlations and experimental data 
exist. This is a matter of validation, which is not the main concern here. Instead, the main concern is 
verification: what are the expected results for a given model, properly implemented, on a refined grid? 
If a model has been implemented consistently in two different CFD codes, then both codes should yield 
the same grid-resolved solution. 

For the flat plate, all results were solved using the same series of structured grids (solved both as 
quadrilateral elements and cut into triangular elements in the unstructured code FUN3D). Grid sizes 
ranged all the way from the finest 545 x 385 to the coarsest 35 x 25, and all were members of the same 
grid family, achieved by removing every other point in each coordinate direction for each coarser grid 
level. The finest grid had a wall normal spacing of y = 5 x 1CT 7 , yielding an average y + of approximately 
0.1 over the plate. The coarsest grid gave average y + « 1.7. The grid extended from x = -0.3333L 
upstream of the start of the plate to x = 2L at the end of the plate. The upper boundary was located at 
y = L. For this and all other reported results, the CFD codes were converged iteratively such that the 
density residual was driven well below I0~ n . 

Fig. 1 (a) and (b) show wall surface skin friction at x = 0.97 and the integrated drag coefficient, 
respectively, using the SA model. Results are plotted as a function of h, a measure of the average 
grid spacing. As h -> 0 (indicating finer and finer grid), results from the two codes and results using 
different grid element types approach each other, demonstrating consistent results as the grid is refined. 
A similar set of plots is shown in Fig. 2(a) and (b) for the SST-V model. 
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Uncertainty estimation procedures are sometimes difficult to use, particularly when (1) CFD results 
exhibit oscillatory convergence behavior, or (2) when results lie outside of the elusive “asymptotic range” 
of grid convergence where Richardson extrapolation is valid. Also, in practice CFD can exhibit appar- 
ently inconsistent order-property behavior. In other words, although one expects a nominally second- 
order code to show p « 2 (where p is the spatial order of accuracy, and the error is proportional to h p ), in 
practice the computed range can vary widely (e.g., from less than 1 to 3 or more). It is important to note 
that as grids are refined, certain characteristics of a scheme such as first-order treatment of boundary 
conditions, which may not affect the overall order-property behavior for coarser grids, can become the 
dominant error source. Also, as discussed in Diskin et al.[1 9], details about the scheme, grid type, grid 
stretching, and body curvature can all have profound effects on the resulting accuracy. 

Using the uncertainty estimation procedure from Celik et al.[20], the finest three grids for the flat plate 
were used to quantify the grid convergence properties. Results are shown in Table 1 . The computed p 
on the finest grid level came out to between 0.78 and 1.98 for the properties looked at, with the lowest 
order tending to occur for the grids cut into triangles. In this table, the fine-grid convergence index 
(defined in Celik et al.[20]) represents a measure of the fine grid solution error, including a factor of 
safety. It is important to note here that particular solution quantities using two different CFD codes are 
not necessarily equal in accuracy, even when using identical grids. Furthermore, the same code on two 
different grids can converge at different rates. Hence it is necessary to perform an independent grid 
resolution study for each code and grid combination, to insure grid-independent consistent solutions. 

A second case focused on the development of the free planar shear layer following the passing of 
two different streams over a thin plate (see Fig. 3). The smaller inner stream had Mach number near 
M = 0.5, whereas the outer larger stream had a Mach number near M = 0.25. This can also be 
considered as a planar co-flowing jet [21]. The Reynolds number was Re = 50,000 based on the full 
inner jet width (L = 1 here). 

For this case only quadrilateral grids have been used to date, ranging from the finest grid with 
327, 680 grid cells to the coarsest grid with 1280 grid cells. The downstream and upper boundaries were 
located at x = 200 and y = 100, respectively. Typical results are shown in Figs. 4 and 5. Again, the 
same turbulence model (SA or SST-V) was shown to approach consistent results for the two codes 
as the grid was refined. Note that some of the cases exhibited oscillatory convergence behavior for 
certain properties. If only one code was being used, this behavior would make it difficult to accurately 
estimate the uncertainty, because an order p cannot be established. For example, Eca and Hoekstra 
[22] assign an uncertainty of 3 times the maximum difference between all available solutions when 
oscillatory convergence occurs. However, since both codes here are consistently approaching the same 
result with grid refinement, the use of two codes in this case can provide a higher confidence estimate 
of the discretization error levels for a given grid and code. This is a strong argument for using more 
than one (verified) code when conducting CFD grid convergence studies. For example, in Fig. 5(a), 
based on the convergence properties of both codes one can confidently assert that the CFL3D fine grid 
solution is well less than 0.5% in error from an infinite-grid solution. However, taken alone the oscillatory 
convergence of CFL3D for this quantity leaves more room for doubt: the method of Eca and Hoekstra 
(using the finest 3 grids) would estimate fine grid solution error to be greater than 4%. 

The final case described here is a 3D bump in a channel. There is no experiment associated with 
this case; it is purely a CFD verification exercise, chosen because the wall bump shape is analytic and 
smooth. There is no flow separation present. Although not discussed here, a 2D verification bump case 
is also provided on the website. The lower wall is a curved viscous-wall bump extending from x = 0 to 
1.5 at the two sides of the computational domain y = 0 and y = -1, but curved to lie further downstream 
at y locations in-between. The maximum bump height (^-direction) is 0.05. The “2D” definition of the 
bump at the y = 0 plane is: 


= 0.05 
= 0 


smf^-^-)] 4 0.3 < x <1.2, y = 0 

\0.9 3.0/ J - - , y 

0 <x < 0.3, 1.2 < x <1.5 


( 11 ) 


The x-location of any position on the bump varies in the spanwise direction between y = 0 and y = - 1 
according to: 


x = xo + 0.3(sin(7n/)) 4 — 1 < y < 0 (12) 

where x 0 is any given ^-location of the “2D” shape at y = 0. A sketch of the 3D bump shape is shown 
in Fig. 6(a). The freestream conditions were M = 0.2 and Re = 3 million per unit length. 
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All results were solved using the same series of structured grids (only quadrilateral grids have been 
used to date). Grid sizes ranged from the finest 65 x 705 x 321 (14.7 million grid points) to the coarsest 
5 x 45 x 21 (4725 grid points). The finest grid had minimum spacing at the wall of y = 1 x 10 -6 , giving 
an approximate average y + of 0.12 over the plate at the Reynolds number run. The coarsest grid had 
y + w 2.0. The grid inflow, outflow, and top boundaries were located at x = -25, x = 26.5, and 2 = 5, 
respectively. 

Only results using SA are shown here. SST-V results using the same grid sequence are similar. 
Pressure coefficients on the bump surface as well as at y = 0 are shown on the 33 x 353 x 161 grid 
in Fig. 6(b) (results are visually essentially identical using both codes on the finer grids). Convergence 
with grid refinement is shown in Fig. 7 for the bump forces, including lift coefficient, drag coefficient, 
and pressure and viscous components of drag coefficient. The two codes produced similar trends and 
approached essentially identical results as the grid was refined, providing further evidence of consistent 
model implementation. 

It should be borne in mind that 3D grid refinement studies that use successive coarsening by remov- 
ing every other point in each coordinate direction very quickly produce grids that are too coarse to yield 
adequate solutions. For example, the 17 x 177 x 81 grid (approximately 244, 000 grid points) produced 
drag levels that were nearly 10% in error from an extrapolated infinite solution, and the 9 x 89 x 41 grid 
(approximately 33, 000 grid points) had error levels near 50%. In order to obtain 3 or more useful and 
usable “fine enough” 3D grids, 3D grid studies are sometimes performed by attempting to parametrically 
employ less than a factor of 2 change in grid spacing in each coordinate direction for successive grids. 
For complex configurations, however, it can be difficult to ensure that the grids produced in this way are 
of the same family. 

4. Conclusions 

A turbulence model verification effort was described, which is part of a larger goal of an AIAA dis- 
cussion group to help improve the consistency, verification, and validation of turbulence models within 
the aerospace community. A NASA website - http://turbmodels.larc.nasa.gov - has been established 
to provide the CFD community with a resource for (1) finding turbulence models, and (2) verifying their 
own coding of turbulence models. This latter capability is made possible by providing simple test cases 
and grids, along with sample results (including grid convergence studies) from one or more previously- 
verified codes. Furthermore, by listing accepted versions of the turbulence models as well as published 
variants, this website establishes naming conventions in order to help avoid confusion when comparing 
results from different codes. 

This paper described some of the turbulence model variations that exist for the popular Spalart- 
Allmaras and Menter shear-stress-transport models, and used two specific versions of these models 
in two independent computer codes to demonstrate consistency with grid refinement for several cases. 
Combining consistency in model implementation in two or more codes along with thorough grid con- 
vergence studies for a range of problems is a powerful methodology that can help to establish the 
grid-resolved solutions and the uncertainties on finite grids. This can be very useful for the common 
situation in which traditional uncertainty methods (such as those involving Richardson extrapolation) fail 
for a single code because of oscillatory convergence or inability to achieve the asymptotic range of grid 
convergence. 
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Table 1. Computed convergence behavior for flat plate case 


Property 

Code 

Computed apparent order, p 

Fine-grid convergence index, % 

Cf at x = 0.97 

CFL3D, SA 

1.98 

0.017 


FUN3D, SA, quad 

1.34 

0.028 


FUN3D, SA, tri 

0.85 

0.592 


CFL3D, SST 

1.21 

0.277 


FUN3D, SST, quad 
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FUN3D, SST, tri 
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CFL3D, SA 
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FUN3D, SA, quad 
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FUN3D, SA, tri 

0.78 
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1.34 
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Figure 1 . Results for flat plate using two independent codes with the SA model; (a) skin friction coefficient at x=0.97, (b) 
drag coefficient. 
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Figure 2. Results for flat plate using two independent codes with the SST-V model; (a) skin friction coefficient at x=0.97, 
(b) drag coefficient. 


7 of 9 


2 - 


1.5 - 


1 - 


0.5 - 


-0.5 - 


Pt/P .=1.046 
Tt/T ref =1.0 

1 quantity from interior 


adiabatic solid wall 


/ 


Pt/P .=1.1862 
Tt/T ref =1 .05 

1 quantity from interior 


y 

adiabatic solid wall 


A 


symmetry 


-10 


Figure 3. Description of planar shear case. 
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Figure 4. Results for planar shear using two independent codes on a quadrilateral grid with the SA model; (a) drag 
coefficient on thin plate, (b) u/a ref along the centerline at x = 95.5. 




Figure 5. Results for planar shear using two independent codes on a quadrilateral grid with the SST-V model; (a) drag 
coefficient on thin plate, (b) u/a ref along the centerline at x = 95.5. 
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Figure 6. Overview of 3D bump test case; (a) sketch of bump shape, (b) pressure coefficients using SA model on 

33 x 353 x 161 grid. 
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Figure 7. Results for 3D bump using two independent codes with the SA model; (a) C L , (b) C D , (c) C D}P , (d) C DjV . 
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